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Abstract 

The popular method of Nose and Hoover to create canonically distributed positions 
and momenta in classical molecular dynamics simulations is generalized to a genuine 
quantum system of infinite dimensionality. We show that for the quantum harmonic 
oscillator, the equations of motion in terms of coherent states can easily be modified 
in an analogous manner to mimic the coupling of the system to a thermal bath 
and create a quantum canonical ensemble. Possible applications to more complex 
systems, especially interacting Fermion systems, are proposed. 
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1 Introduction and summary 



The typical problem in statistical physics is the determination of ensemble av- 
erages. The canonical ensemble is characterized by a constant temperature, i. 
e. the total energy of the system is allowed to fluctuate around its mean value, 
but the system is kept at a constant temperature by thermal contact with an 
external heat bath. Besides the direct evaluation of ensemble averages which is 
impossible in many cases, especially in interacting many-body systems, numer- 
ous different approaches have been developed to calculate canonical ensemble 
properties, among them Monte Carlo approaches and stochastic techniques. 
In classical molecular dynamics, Nose has developed a scheme that allows 
to calculate canonical averages by averaging over a deterministic isothermal 
time evolution [1,2]. This scheme is called the classical Nose-thermostat and 
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has attracted much interest. In the Nose-method, a degree of freedom s and 
its conjugate momentum p s are added to the original system for temperature 
control. The additional degree of freedom s acts as a scaling factor for the 
positions and momenta of the original system. The idea is now that the isoen- 
ergetic time evolution of the extended system (that conserves the total energy 
of the extended system) yields an isothermal time evolution in the subspace 
of the variables of the original system. This holds for ergodic time evolutions. 
For a detailed review, see [2]. 

Although in practice the original formulation of Nose turned out to be too 
cumbersome and featured ergodicity problems in many cases, it allowed for 
a number of improvements that led to very effective and versatile methods 
[3]. Simply speaking, the resulting schemes exploit the equipartition theorem 
of classical mechanics to determine the equations of motion of pseudofriction 
coefficients. The most reliable methods are the so-called Nose-Hoover chains 
[4] and the demon method of Kusnezov, Bulgac, and Bauer (KBB) [5] for 
which even the simple one-dimensional harmonic oscillator is ergodic. 

For quantum systems, equivalent methods of comparable power are not yet 
available. Grilli and Tosatti have found a theorem that provides a basis for 
a seemingly possible translation of the Nose-method to quantum mechanics 
[6]. However, in practice, their method features substantial problems [7,8]. 
Kusnezov has proposed a method for finite dimensional quantum systems 
that can be applied if all eigenvectors and eigenvalues of the Hamiltonian are 
known [9]. 

For coherent states, a quantum phase space (r, p) can be defined properly and 
a thermal weight function w qm ((3;r,p) exists that permits the calculation of 
canonical ensemble averages as phase space integrals [10]. In this article, we 
present a modification of the quantum equations of motion of coherent states 
in a one-dimensional harmonic oscillator, following closely the ideas of Nose, 
Hoover, and KBB for classical systems. In order to calculate ensemble averages 
by time averaging, the quantum equations of motion of the parameters (r,p) 
of the coherent states are modified. More precisely, a classical pseudofriction 
coefficient p v is added to the system and the equations of motion are designed 
in such a way that the distribution 



which is defined on the mixed quantum-classical phase space (r,p,p v ) is a 
stationary solution of a generalized Liouville equation. As a consequence, if 
the system is ergodic, / is the stationary probability distribution generated 
by the modified quantum dynamics, and canonical ensemble averages can be 
calculated by time averages over the trajectories thus generated. Hence, our 
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method provides an isothermal quantum time evolution for a quantum system 
of infinite dimensionality. It is straightforward to generalize it to systems of 
many distinguishable particles in a three-dimensional harmonic oscillator as 
well as to free particles. Even non-interacting fermions may be thermalized, 
since in this case, the quantum distribution function is also known [11]. 



2 Method and setup 

2. 1 Coherent states in a harmonic oscillator potential 

Given the Hamilton operator H of the one-dimensional harmonic oscillator 

H = hu (aW^) » (!) 
coherent states are defined as eigenstates of the destruction operator a 
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A coherent state is labelled by its complex eigenvalue z which corresponds to a 
pair of real parameters (r,p). Explicitly, in coordinate representation, coherent 
states are shifted Gaussian wave packets characterized by the parameters r 
(mean position) and p (mean momentum): 



(2) 



. . . . f (x — r) 2 mw i \ . . 

{x\z) = {x\r,p} ocexp<^ ^- + ^px'>. (3) 

Coherent states have been extensively investigated [12]. In particular, the fol- 
lowing equality is useful for considering the time evolution of coherent states 
in a harmonic oscillator potential: 



exp(— iuajat) \z) — \ exp(— iut)z ) . (4) 

This implies that the exact quantum time evolution of a coherent state in a 
harmonic oscillator potential is given by the following equations of motion for 
the parameters r and p 



We stress that in a harmonic oscillator potential the solution of these two cou- 
pled ordinary differential equations provides the exact quantum time evolution 
of coherent states. 

Furthermore, the set of all coherent states forms an overcomplete basis of 
the Hilbert space with / j0£ \r,p)(r,p\ = 1 As a consequence, given an 
observable B, its thermodynamic mean value may be evaluated using coherent 
states: 



1 f dr dp I rj —/3H | v 

= ^)i(M) (r ' p| ^ e Hr ' p) ' 

where (3 = is the inverse temperature, Z((3) = tr ( e_/3 ~) is the usual 
canonical partition function and ((■)) denotes canonical averages. 

As shown in [10], one can interpret the space of the continuous parameters 
r and p as a phase space and rewrite (6) as a phase space integral with the 
thermal weight function w qm (P;r,p) 
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dr dp 
Z{(3) J (2nfi) 
' drdp 



w qm {(3]r,p) B(r,p) , 



Wqm(/3;r,p) , 



(7) 
(8) 



with 



B(r,p) = (r,p \ B \ r,p) , (9) 
w qm (P;r,p) = e-^ 2 (^-i) = e -£+W(^-i)/(fc) . (10) 

The function w qm contains all quantum statistical properties of the system. 
From (10) it can be inferred that formally, it differs from the classical dis- 
tribution function of the harmonic oscillator by the factor [eP^ — 1^ jifihS). 
Note that this factor tends to 1 in both the classical (h~ — > 0) and the high- 
temperature (P — > 0) limit. 



2.2 Modification of the equations of motion 



The idea of our method is to modify the equations of motion (5) of the coherent 
states in such a way that the distribution function w qm is sampled provided 
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the time evolution is ergodic. To this end, we proceed in a way which is in close 
analogy to the approaches in classical molecular dynamics [3-5]. The equation 
of motion of the parameter p is supplemented by a term similar to a frictional 
force. The time evolution of the pseudofriction coefficient is then determined 
by the condition that the desired distribution function is a stationary solution 
of a generalized Liouville equation in the generalized phase space. 



2.2.1 Nose-Hoover thermostat and Nose-Hoover chain 

Adopting the notation of Martyna et al. [4], we investigate the following ana- 
logue of the classical Nose-Hoover dynamics for the quantum dynamics of 
coherent states: 



d p d 2 Pv m\ 

—r = — , — p = —muo r — p—r . (11) 
dt m ' dr Q y ' 

The key point is the equation of motion of the pseudofriction coefficient p v . It 
is determined by the condition that the distribution function 



p 2 

f((3;r,p,p v )(xw qm ((3;r,p) ex P(-^^g) ( 12 ) 

/ p 2 1 2 2 f^-l p 2 v \ 

oc exp — ( — muo r ) — P—h 

l \ v 2m 2 ' frw 1 2Q) 

is a stationary solution of the following generalized Liouville equation in the 
mixed quantum-classical phase space Y = (r,p,p v ): 



= -/ • (— r + —P+ —p ) 

\dr T dp^ dp v ^ v J 

We calculate the left-hand side of (13), employing the equations of motion of 
r and p, (11): 



dt dp^ dr dprj^ 



(14) 
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On the right-hand side of (13), we have the freedom to impose the constraint 
dp v /dp v = that is common in this context [2]. We obtain 




-f%. (15) 
Equating (14) and (15) yields the following equation of motion for p v : 



d 1 (p 2 f» - 1 \ 

Again, the only difference between this equation and its classical counterpart 
is given by the factor [eP^ — I) / {(3?kij) . Moreover, (16) retains the property 
of its classical counterpart that the time evolution of the pseudofriction coef- 
ficient is governed by the deviation of the actual value of a quantity related 
to the kinetic energy from its canonical average value. This can be inferred by 
evaluating 



2 phw _i 
m noj 11 
using (7). 

Finally, it is easily confirmed that the set of dynamical equations (11), (16) 
conserve the quantity 




The equations (11) and (16) form a genuine quantum Nose- Hoover thermostat 
for coherent states. Since in classical molecular dynamics, these equations of 
motion frequently feature ergodicity problems, Martyna et al. have developed 
the idea of a chain thermostat [4] . This method implies to impose another ther- 
mostating pseudofriction coefficient on p v which may be coupled to yet another 
pseudofriction coefficient, and so on, thereby forming a chain of thermostats. 
The application of the idea to the quantum case does not infer anything new 
compared to the classical case, since only the first pseudofriction coefficient 
of the chain interacts with the quantum phase space variables. Therefore, for 
further particulars we refer the reader to [4]. 
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2.2.2 KBB-method 



Another generalization of the Nose-Hoover thermostat that is frequently used 
in classical molecular dynamics is the so-called demon method proposed by 
Kusnezov, Bulgac, and Bauer [5]. The advantage of this method is that the 
Hamilton function of the envisaged system does not have to contain a kinetic 
energy term for temperature control; instead, the time derivative of the tem- 
perature control variables is postulated to be proportional to the difference of 
two arbitrary quantities whose ratio of canonical averages is 1//3. 

At least two pseudofriction coefficients, so-called demons, are introduced for 
temperature control. Both the equations of motion for positions and momenta 
are supplemented by additional terms. We introduce the demons into the 
quantum equations of motion of the parameters of coherent states: 



^r= ^-g' 2 (QF(r,p) , j t p = -muj 2 r - g[(()G(r,p) . (19) 

F(r,p),G(r,p) are arbitrary functions of the quantum phase space variables. 
<7i(C)>02(£) are functions of the demons which have to be chosen so that the 
integration of the distribution function / converges. g[,g 2 are the respective 
derivatives. The distribution function on the phase space T = (r,p,(,0 reads 



f(r,p,£,0 = exp - — + -mu r — (3{ + 20 

\ 2m 2 nw k 2 «i / 

and the time evolution of the demons is, as above, deduced from the require- 
ment that / is a solution of a generalized Liouville equation in the phase space. 
We obtain 



d ( p e f»>"-l ldG\ 

d , (2 r.e^-l ldF\ 

Again, it is interesting to notice that 



i.e. the ratio of the canonical averages of the quantities that determine the 
time derivative of the demons is j3, just as in the classical case. The quantity 
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p 2 1 2 2 , e^-l <? 2 (0 <?i(C) , 0A s 

H = — muj r ) — — 24 

v 2m 2 ; ki /t 2 v ; 



is conserved during the time evolution defined by (19), (21), (22). 

In principle, since the choice of the functions F,G,gi,g 2 is arbitrary, this 
method offers a lot of freedom. The most prominent coupling scheme recom- 
mended by KBB is the so-called cubic coupling scheme with the following 
choice of functions [5] 



9i = \i 2 , 92 = \C , F = r 3 , G = p , (25) 
which leads in the quantum case to the special set of equations of motion 



d p 

—r = £r 3 , — p = —mto 2 r — ( d p (26) 

at m 




l^-K^-H <28) 

that we have investigated taking k± — k 2 — 1. Finally, we note that (27), 
(28) may easily be linked to the equations of motion proposed by Kusnezov 
in [9]. w qm plays the role of Kusnezov's p(Q,P). However, while Kusnezov's 
approach is limited to quantum systems of finite dimensionality, our method 
works for this system with a Hilbert space of infinite dimensionality because 
we take advantage of the properties of coherent states. 



3 Results 



In classical molecular dynamics simulations of the harmonic oscillator, the sim- 
ple Nose-Hoover method features ergodicity problems, while the Nose-Hoover 
chain method and the demon approach of KBB work well [4,5]. The correct 
classical phase space density is perfectly reproduced both by the chain and by 
the demon dynamics. 

Formally, the only difference between the quantum phase space density w qm , 
eq. (10), and its classical counterpart is given by the factor [eP^ — 1) /((3hu). 
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This also applies to the respective equations of motion of the different dynam- 
ics. Since this factor is only a number that depends on temperature, but not 
on the phase space variables r and p, we anticipate that it does not influence 
the overall characteristics of the dynamics. Therefore, we expect that ergod- 
icity problems in the quantum case will arise under the same circumstances 
as in the classical case. 

We show results for the set of parameters chosen in [4] to enable a direct 
comparison. We took m — \,u — 1 and initial conditions r(0) = l,p(0) = 1 
with 1/(3 = 1.0. The numerical integration of the equations of motion was 
carried out with a fourth-order Runge-Kutte algorithm with a step size that 
ensured conservation of pseudoenergy to more than seven significant figures. 
All runs were made over a total integration time of 2000r, where r = 1-k/uj. 



3.1 Nose-Hoover and Nose-Hoover chain method, 
KBB method with cubic coupling scheme 



The left panel of figure 1 presents a (r, p)-density map and the projected dis- 
tribution functions for the simple Nose-Hoover dynamics. We find a result that 
is very similar to the classical case [4]: The dynamics does not fill the phase 
space with the correct weight and, moreover, we find that the obtained distri- 
butions strongly depend on the initial conditions and values of the parameters 
chosen (not shown). Thus, the dynamics is not ergodic. 

The situation changes radically with the introduction of a second thermostat- 
ing variable acting upon the first pseudofriction coefficient p v , see figure 1, 
right panel. The distribution functions sampled by this time evolution repro- 
duce w qm extremely well, and changes of the initial conditions and parameters 
do not have a noticeable effect on the results. The dynamics generated in this 
way is obviously ergodic. The addition of further thermostating variables does 
not influence the results. 

We point out that the statistics obtained by time averaging over the quantum 
time evolution is the quantum statistics of the harmonic oscillator. To make 
this evident, we present plots of the partially integrated distribution function 
w qm (r,p)/Z along with plots of its classical limit 



1 / p 2 1 \ 

w d (r, p) = Hm — — w qm (r, p) = (3hw exp -/?(— + -mw 2 r 2 ) (29) 
e^-i_! Z{(3) \ 2m 2 J 

which is proportional to the classical canonical ensemble distribution function. 
Since {eP^ — / ((3Huj) > 1 for all (3, the quantum distribution function is 
always narrower compared to its classical limit. 
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Fig. 1. Left panel: Simple Nose-Hoover dynamics of a quantum harmonic oscil- 
lator, right panel: Nose-Hoover chain dynamics. From above: (r, p)-density plot, 
momentum distribution, position distribution. The solid line depicts the exact 
quantum result given by the respective partially integrated function w qm /Z (e.g., 
f(r) = \ J j£ h w qm (r, p) ) , the dashed line represents the corresponding classical 
distribution w c i normalized to the same value (see (29)). The distributions sampled 
by time averaging are presented as histograms. 

Figure 2 presents the results obtained from a KBB-demon-dynamics using 
the cubic coupling scheme. The results are similar to the case of the chain 
dynamics, in particular, the dynamics is also ergodic. 

3.2 Mean values of selected observables 

Finally, the results of time averaging are compared to the analytical ensemble 
averages for two typical observables, the internal energy and its variance. The 
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Fig. 2. KBB dynamics of a quantum harmonic oscillator, description as in figure 1. 
analytical formulas are briefly given: 

^ = «5-»-«5» , = ( 5S 5^)' 

Figure 3 shows the excellent agreement between the exact results and the 
results obtained by time averaging with a KBB dynamics. The small deviations 
are clearly of statistical origin. They increase at higher temperatures because 
we kept the total sampling time constant, although the volume of the relevant 
phase space increases with temperature. Therefore, to cover it with the same 
accuracy at a higher temperature, one would need a longer sampling time. 
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Fig. 3. Values of the internal energy and its variance of the harmonic oscillator 
obtained from time averaging with a KBB dynamics (crosses) compared to the 
exact quantum canonical ensemble result (solid line). 

4 Discussion and outlook 

This article presents a straightforward, yet non-trivial extension of the power- 
ful methods of heat bath coupling in classical molecular dynamics simulations 
to a genuine quantum system of infinite dimensionality. The application of 
the method to a quantum system of many distinguishable particles or a three- 
dimensional harmonic potential is a simple generalization. 

Since the knowledge of w qm is indispensable for the setting up of the equations 
of motion for the pseudofriction coefficients, the method is limited to systems 
where w qm is known. Therefore, also non-interacting identical fermions, both 
moving freely or contained in a harmonic oscillator potential, can be thermal- 
ized using the respective distribution functions [11]. Moreover, by coupling 
one of the solvable systems to a more complex system of interacting particles, 
one can possibly determine its equilibrium properties. This idea that permits 
to evaluate ensemble averages by time averaging is potentially very power- 
ful, since efficient approximate quantum dynamics methods (Time-Dependent 
Hartree-Fock, Fermionic Molecular Dynamics, etc. [13]) are available that are 
applicable also for indistinguishable fermions. Thus, a new method of calcu- 
lating thermodynamic properties of interacting Fermion systems seems con- 
ceivable that might also work where other methods fail. 
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